rm(list=ls())

# Set working directory 
setwd("")

# Initialize stuff
rlnm <- c("Log", "Brier", "Spherical", "Boosting", "As1", "As2")
preds <- replicate(6, matrix(0, 100, 3), simplify = F)
csf <- replicate(4, matrix(0, 99, 6), simplify = F)
c <- seq(from = 0.01, to = 0.99, by = 0.01)
x <- 2

# Loop over DGPs
for (dgp.nr in 1:4){
  
  if (dgp.nr <= 2){
    xlim <- c(-2.5, 2.5) 
  } else if (dgp.nr == 3){
    xlim <- c(-1, 4)
  } else if (dgp.nr == 4){
    xlim <- c(0, 10)
  }
  xseq <- seq(from = xlim[1], to = xlim[2], length.out = 40)

  # Load workspace -> NOTE: These files are produced by running "MCstudy_rollingwindows.R"
  load(paste0("rwMC_dgp", dgp.nr, "_120_10000.RData"))
  
  # Unlist theta
  theta.all <- sapply(1:6, function(j) sapply(pars, function(s) s[, j]))
  pdat <- matrix(0, length(xseq), 7)
  pdat[,1] <- dat$p.fct(xseq)
  for (j in 1:6){
    pdat[, j+1] <- rowMeans(outer(xseq, theta.all[, j], function(x, y) logit(x*y)))
  }
  
  # Prediction plot (mean over MC iterations)
  mccolors <- colorRampPalette(c("blue","red","green"))
  pdf(paste0("pred", dgp.nr, "_rw.pdf"), family = "Bookman")
  par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
  matplot(x = xseq, y = pdat, type = "l", col = c(1, mccolors(6)), lty = c(rep(1,4), rep(2,3)), lwd = 4, 
        xlab = "x", ylab = "")
  legend("topright", inset = c(-0.4, 0), c("True", "Log", "Brier", "Spherical", "Boosting", "As1", "As2"),
       col = c(1, mccolors(6)), lty = c(rep(1,4), rep(2,3)), lwd = 4, xpd = TRUE)
  dev.off()
  
  # Classification plot (mean over MC iterations)
  pdf.options(family = "Bookman")
  pdf(paste0("class", dgp.nr, "_rw.pdf"), family = "Bookman")
  tr <- (dat$p.fct(x) > c)
  gd <- t(sapply(c, function(l) colMeans(logit(theta.all*x) > l)))
  par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
  matplot(cbind(tr, gd),type = "l", lwd=4, x=c, col=c(1,mccolors(6)), ylab="", xlab=expression(c), 
        lty=c(rep(1,4),rep(2,3)))
  legend("topright", inset = c(-0.4, 0), c("True", "Log", "Brier", "Spherical", "Boosting", "As1", "As2"),
       col = c(1, mccolors(6)), lty = c(rep(1,4), rep(2,3)), lwd = 4, xpd = TRUE)
  dev.off()
  
  # Clear most of the workspace
  rm(list = setdiff(ls(), c("dgp.nr", "theta.all", "preds", "rlnm", "xlim", "xseq", "c", "x", "logit", "csf", "dat")))
  
}